// Simple CUDA example by Ingemar Ragnemalm 2009. Simplest possible?
// Assigns every element in an array with its index.

// nvcc simple.cu -L /usr/local/cuda/lib -lcudart -o simple

#include <stdio.h>

const int N = 1024; 
const int size_x = 512;
const int size_y = 1;
const int X = 64; 
const int Y = 32; 

__global__ 
void add_matrix(float *a, float *b, float *c, int N)
{
  int x = threadIdx.x + blockIdx.x * blockDim.x;
  int y = threadIdx.y + blockIdx.y * blockDim.y;
  int index = x+y*blockDim.x*gridDim.x;
  //  int index_no_coal =  y+x*blockDim.y*gridDim.y;

  c[index] = a[index] + b[index];
  //c[index_no_coal] = a[index_no_coal] + b[index_no_coal];
}


int main(int argc, char **argv)
{
  //N = atoi(argv[1]);
  dim3 dimBlock( size_x, size_y );
  dim3 dimGrid( X, Y );

  //float a[N*N];
  //float b[N*N];
  //float c[N*N];

  float *a = (float *) malloc(sizeof(float) * N * N);
  float *b = (float *) malloc(sizeof(float) * N * N);
  float *c = (float *) malloc(sizeof(float) * N * N);

  float *a_g;
  float *b_g;
  float *c_g;
  const int size = N*N*sizeof(float);

  float theTime;
  cudaEvent_t startEvent,stopEvent;
  cudaEventCreate(&startEvent);
  cudaEventCreate(&stopEvent);

  
  cudaMalloc( (void**)&a_g, size);
  cudaMalloc( (void**)&b_g, size);
  cudaMalloc( (void**)&c_g, size);
  
  
  for (int i = 0; i < N; i++){
    for (int j = 0; j < N; j++){
      a[i+j*N] = 10 + i;
      b[i+j*N] = (float)j / N;
      c[i+j*N] = 0;
    }
  }
 
  cudaMemcpy( a_g, a, size, cudaMemcpyHostToDevice );
  cudaMemcpy( b_g, b, size, cudaMemcpyHostToDevice ); 
  cudaMemcpy( c_g, c, size, cudaMemcpyHostToDevice ); 

  cudaEventRecord(startEvent, 0);
  add_matrix<<<dimGrid, dimBlock>>>(a_g,b_g,c_g,N);  
  cudaEventRecord(stopEvent, 0);
  cudaEventSynchronize(stopEvent);  
  cudaEventElapsedTime(&theTime, startEvent,stopEvent);

  cudaMemcpy( c, c_g, size, cudaMemcpyDeviceToHost );
  cudaFree(a_g);
  cudaFree(b_g);
  cudaFree(c_g);
  
  for (int i = 0; i < N; i++)
    {
      for (int j = 0; j < N; j++)
	{
	  printf("%0.2f ", c[i+j*N]);
	}
      printf("\n");
    }
  printf("time: %0.3f ms\n", theTime);
  
  free(a);
  free(b);
  free(c);
  return EXIT_SUCCESS;
}
